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1. Introduction 

The Higgs boson is an elusive particle which couples to particles according to their mass and 
so is weakly coupled to quarks and leptons. The dominant production mechanisms of Higgs 
bosons at hadron colliders is from gluon-gluon fusion and vector boson fusion. In this paper, 
we consider another mechanism which may be more relevant to an experimental search. 
This is the associated production of Higgs bosons with vector bosons. This production 
mechanism, pp — > WH/ZH + X, is the most promising discovery channel for a light 
Standard Model (SM) Higgs boson at the Tevatron. This is because the Higgs, which 
decays predominantly into bb pairs, can be tagged by the associated vector boson. 

Such processes can be simulated in parton shower generators which resum soft and 
collinear leading logarithmic, as well as an important subset of next-to-leading logarithmic 
contributions to all orders. These simulations can further be improved by matching the 
parton shower to higher order matrix elements. One way in which this is done in the generic 
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Figure 1: Associated WH boson production. 



parton shower generators is through the use of the matrix element correction [1] which 
generates harder emissions in regions outside the reach of the parton shower at a rate given 
by the matrix element. A more rigorous matching procedure is the MCONLO method [2-6] 
which has been implemented for a multitude of processes in the HERWIG event generator [7] 
and for some processes [8] in its successor the Herwig++ event generator [9,10]. More 
recently, another matching method was proposed, called the POsitive Weighted Hardest 
Emission Generator (POWHEG) [11,12], which achieves the same aim as MCONLO, with the 
creation of positive weighted events and is furthermore independent of the shower generator 
used. The POWHEG method has been applied to Z pair hadroproduction [13], heavy flavour 
production [14], e + e~ annihilation to hadrons [15], Drcll-Yan vector boson production 
[16, 17] and top pair production at the ILC [18]. 

In this paper, we aim to simulate the NLO hadronic decay of the light Higgs boson 
produced in association with a vector boson using the MC@NL0 method. The parton shower 
generator we will be employing is Herwig++. In Section 2, we first discuss the cross- 
sections for associated production of the Higgs boson with a vector boson and its subsequent 
hadronic decay rate. We then discuss the application of the MC@NL0 method to the decay in 
Section 3. In Section 4, we show some comparative distributions obtained from the parton 
shower and in Section 5 we summarize our conclusions. Finally, it should be noted that in 
this paper, we do not apply the MC0NL0 method to the initial state emissions. 

2. Cross-sections and decay rates 

2.1 Associated Higgs production with a W boson from qq annihilation 

The process q(p q ) + q' (p'q) — ► W* — ► W(pw) + H(pn) is illustrated in Figure 1. If we define 
the center of mass energy squared of the partonic system by 



s = 



(Pq + Pq'f , 



(2.1) 



we have for the differential cross-section [19], 



dcosO* 




cos 2 6 w 48vrs 3 /2 ^ s - 




[2(1 -/^)+/^ sin 2 0*], (2.2) 
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Figure 2: Associated ZH production. 



where 9* is the angle between the W boson and the quark in the partonic CMF, Gf is the 
Fermi coupling constant, &w is the Weinberg mixing angle, Mw and Mu are respectively 
the W and Higgs boson masses and V q qi is the relevant CKM matrix element. (3w is the 
speed of the W boson in the partonic CMF and is given by 

y/[s - (M w + M H y\[s - (M w - M H Y] 



s-M 2 H + M w 



(2.3) 



The relativistic boost factor is (1 — 0^) 1 I 2 . Equation 2.2 when integrated gives for 
the total partonic cross-section, 



GlM^V 2 -, (3 w7w M w (** 
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2/3; 



(2.4) 



Convolving this with parton distribution functions (PDFs), we obtain the hadronic cross- 
section as 



= dxi / dx 2 [f g (x 1 ,Mw)fq>(x2,Mw) + x 1 -^ x 2 ](Tp, (2.5) 

Jt Jt/xi 

where x±,x 2 are the momentum fractions of the incoming partons and taking S as the 
hadronic beam-beam center-of-mass energy, we have s = x±x 2 S. 

2.2 Associated Higgs production with a Z boson from e + e annihilation 

The differential cross-section for the process e + e~ — > Z* — > ZH is given by [20] 
da _ GpMyy PzizM* 



dcos9* 



M 2 



Ml 



Mf 



(2.6) 



COS 2 #VK 32vrs 3 / 2 

[1-4 sin 2 9 w + 8 sin 4 6 W ] [2(1 - /?§) + /?| sin 2 0*] , 

where 9* is the angle between the electron and the Z boson in the CMF, s is the center-of- 
mass energy and 7^ and (5z are obtained from the analogous expressions for W production 
in section 2.1 by substituting Mz for Mw- Equation 2.6 integrated over cos 6* gives 

GlM 2 &zlzM\ fs + M 2 z ' 2 
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Figure 3: Lowest order Higgs hadronic decay rate. 
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Figure 4: Virtual corrections (a): vertex correction (b),(c): self-energy corrections. 

2.3 Higgs boson decay to bb pairs 
2.3.1 Lowest Order decay rate 

The lowest-order decay rate for this process, illustrated in Figure 3, is given by 



Tb(H — > bb) 



3GFml -M H pi, 



(2.8) 



where nib is the mass of the bottom quark and (3q = \/l — -^r- 



2.3.2 Virtual radiative corrections 



If one uses the on-shell renormalization scheme, the self-energy diagrams in Figures 4(b) 
and 4(c) are cancelled by counter-term diagrams leaving us with the vertex correction in 
Figure 4(a) and its counter-term. When evaluated in the massive gluon regularization 
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Figure 5: Real gluon emission. 



scheme, the final result is [21], 
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(2.9) 



where we have introduced the gluon mass fi to regulate the infrared singularity in Figure 
4(a). 



2.3.3 Real emission corrections 



Using the following relations for the energy fractions of the b and b quarks in terms of the 
parton momenta in the Higgs rest frame, 



Xb 



2Ph ■ k 
Ml ' 

2 Pb -k 
Ml ' 



(2.10) 



we have for the real emission decay rate, the following expression: 
dT R _, a s C F 



dx b dx b 



= IV 
= r B 
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where p 



Mi 



This can be integrated in the massive gluon scheme to get, 
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Summing this with Tb and IV in equation 2.9, the dependence of the gluon mass p, disap- 
pears to give, 



where 
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Also from equations 2.9, 2.11 and 2.12, we note that Ty can be written in terms of the real 
emission matrix element squared M as 
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3. MC@NL0 method 

The phase space for gluon emission is given in terms of the Dalitz plot variables Xb,x b by 

(3.1) 



\(xl - p, x\ - p, {2-x b - xif) < , 
where the function A is defined by, 

X(x, y, z) = x 2 + y 2 + z 2 — 2xy — 2yz — 2xz . 



(3.2) 
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Figure 6: Phase space for gluon emission. 



This is equivalent to the condition 

(1 - x b )(l - xi)(x b + xi - 1) > p(2 -x b - x b ) 2 , 



(3.3) 



Figure 6 shows the corresponding phase space region where we have labeled as P the 
emission region and O the region outside the P but in the half-triangle (1 — x b )(l — 
x b ){x b + x b - 1) > 0. 

Now using equation 2.11 and 2.15, we can re- write T^lo i n integral form as, 



asCF {M-AU v } + ^M 



2vr 



2vr 



+T B I dx b dxi 
o 



rNLO =<TB dx b dx b 

(3.4)" 

Now, if we define a functional Ti which represents hadronic final states generated by a 
parton shower starting from a configuration i, we can write down an overall generating 
functional for hadrons from Higgs boson decay as 



dx b dx b 



+ T B / dx b dx b F bb 
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(3.5) 



where T bb is the functional representing shower final states resulting from the process 
H — > 66 and jF fe ^ represents final states from H — > 6651. 

There are two problems with this functional as it is written above. The first is the 
highly inefficient sampling that will be required to generate starting 66 and bbg configu- 
rations according to the A4 since it is divergent in the soft and collinear regions of phase 
space. The second problem arises because when interfaced with the parton shower, leading 
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order configurations starting with bb would radiate quasi-collinear gluons with a distribu- 
tion given by the parton shower approximation to Ai which we shall call Mc- These are 
already included in the starting bbg configurations. Likewise, some of the bbg configurations 
would include 66-like configurations if the gluon is quasi-collinear to either the quark or 
anti-quark. This problem is often referred to as double- counting. 

Before we discuss how to solve these problems, let us investigate as ^ F M-c- In the 
parton shower generator Herwig++, this is the massive quasi-collinear splitting function 
for the emission of a gluon from a quark [22] : 



2vr 



Mcdq 2 dz 



cxsCf dq dz 
2vr q 2 1-z 



1 + z' 



2m z b 
zq 2 



(3.6) 



Here z and q are Herwig++ evolution variables and are respectively the light-cone momen- 
tum fraction of the quark after the emission of the gluon and an angular variable related to 
the relative transverse momentum of the quark after the emission [23]. The splitting func- 
tion in equation 3.6 for the b quark can be re-written in terms of the Dalitz plot variables 
as 

dxbdxi 



2vr 



where if r is defined as 



then z is given by 
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Note that Xb, x b are given in terms of z, r and q by 

^2 



Xb 



1 - z(l - z) 
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= (2 - x b )r + (z- r)^x 2 -Ap . 



(3.7) 



(3.8) 



(3.9) 



(3.10) 



For emission from b anti-quark interchange Xb and x b in the equations above. 

In Herwig++, the quasi-collinear region of phase space covered by the parton shower 
is defined by imposing the condition, 



M 



(3.11) 



The corresponding regions are shown in Figure 7 labeled Jb, Jg whilst the unpopulated 
dead region is labeled D. Note that region P in Figure 6 corresponds to the union of 
regions Jb, Jg an d D. Now going back to the functional defined in equation 3.5, if in the 
region J = Jb U Jg, we subtract the parton shower approximation M.q from the second 
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Figure 7: Phase space showing hard D and soft/collinear Jb, Js ghion emission regions. 



term in equation 3.5 and add it to the first term we get 



T = T B jf dx b dx- b T hl 1 2 - ^^{M -M c - 4n v )| + T x 
+ T B / dx b dx 
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This is the MC@NL0 method which solves the problem of double-counting in the parton 
shower regions J. It also solves the problem of the sampling inefficiency since M — > Mc m 
the divergent regions and therefore M — Mc tends to there. In Appendix A, we describe 
the algorithm used for the evaluation of the above integrals and the generation of events. 
We also discuss how we regularize some residual divergences by the use of mappings in 
Appendix B. 

The procedure followed for event generation for associated Higgs production is outlined 
below. 

1. For pp annihilation use equation 2.5 to distribute the Mandelstam variables x±,X2 
and the angle 9* according to the differential cross-section. From these variables, re- 
construct the W and Higgs boson four-momenta. For e + e~ annihilation, use equation 
2.6 to distribute the angle 9* and reconstruct the Z and Higgs boson four-momenta. 

2. In the rest frame of the Higgs boson, generate the Dalitz plot variables Xb,xi and 
event weight as described in Appendix A. Reconstruct the four-momenta of the 
quark, anti-quark and gluon in this frame. 
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3. Boost the four-momenta back to the lab frame. 
More details about the MCONLO method can be found in [2]. 

4. Results 

Following the prescription above, events for associated Higgs boson production of mass 
114 GeV and decaying into bb pairs with m& = 5 GeV, were generated and interfaced with 
Herwig++ 2.3.0 [10]. The following approximations were considered: 

1. The Herwig++ parton shower interfaced to leading order events (LO), 

2. the parton shower interfaced to leading order events which are supplemented by 
matrix element corrected events in the dead region (ME), 

3. the parton shower interfaced to events generated by the MC@NL0 method (MC@NLO). 
The two processes considered were: 

1. associated WH production from qq' annihilation at the Tevatron (1.96 TeV), 

2. associated ZH production from e + e~ annihilation at the ILC (0.5 TeV). 

The following distributions were considered and shown in Figures 8-10 for qq' annihilation 
and associated W boson production and Figures 11 - 13 for e + e _ annihilation and asso- 
ciated Z boson production. In the simulations, only the leptonic decays of vector bosons 
were considered. 

1. The mass of the bb pair before hadronization, 

2. the energy of the bb pair, 

3. the transverse momentum of the bb pair with respect to the beam axis, 

4. the transverse momentum of the bb pair with respect to the direction of the vector 
boson, 

5. the longitudinal momentum of the bb pair with respect to the beam axis, 

6. the rapidity of the bb pair. 

From the mass and energy reconstruction plots in Figures 8 and 11, we see that the extra 
gluon radiation simulated by the MC@NL0 and matrix element correction methods smooth 
out the respective peaks due to the production of more low mass and energy bb pairs. We 
also observe as expected that the matrix element correction method underestimates the 
amount of hard gluon radiation when compared to the MC@NL0 distributions. 

For qq' annihilation, the effect of this extra radiation on the transverse momenta plots 
is diminished by the initial state radiation from the incoming partons as can be seen in 
Figure 9. The effect is more clearly seen in the corresponding plots for e + e~ annihilation 
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Mass of b anti — b pair (Mbb) Energy of b anti — b pair (Ebb) 




Figure 8: Mass and energy of the bb pair (qq' annihilation). 



pT of b anti — b pair wrt to beam axis pT of b anti-b pair wrt to the W boson 




Figure 9: of the bb pair w.r.t the beam axis and w.r.t the W boson (qq' annihilation). 

in Figure 12 where there is no initial state radiation. At leading order, the bb pairs are 
produced predominantly at right angles to the beam axis (equation 2.6), hence the peak 
at high pt m the LO plot for e + e~ annihilation. The effect of extra gluon emission is to 
smooth out the peak as seen in the MCONLO and ME correction plots. Likewise at leading 
order, the bb pairs are produced back-to-back with the associated Z boson. The effect of 
gluon radiation is therefore to increase the transverse momentum of the bb with respect 
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Figure 10: Longitudinal momentum and rapidity of the bb pair (qq' annihilation). 
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Figure 11: Mass and energy of the bb pair (e + e annihilation). 



to the Z boson. This is what is observed in the MCSNLO and ME correction distributions 
with the latter underestimating the amount of radiation in the tails in comparison with 
the former. 

Finally, in Figures 10 and 13, the MCONLO and ME correction method predict slightly 
more peaked distributions for the longitudinal momenta around the central value of GeV. 
This is a result of the loss of energy from the bb pairs due to extra gluon radiation. The 
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pr of b anti — b pair wrt to beam axis 



pT of b anti — b pair wrt to the Z boson 




Figure 12: px of the bb pair w.r.t the beam axis and w.r.t the Z boson (e + e annihilation). 



Lorigtitudinal momentum of b anti— b pair (Lbt>) 



Rapidity of b anti — b pair (77 bb) 




Figure 13: Longitudinal momentum and rapidity of the bb pair (e + e annihilation). 



rapidity distributions show the MC@NL0 and ME correction methods predict the production 
of more high rapidity bb pairs in comparison to the leading order prediction. This arises as 
a result of more bb pairs being produced at low pp with respect to the beam axis (Figure 
12) and therefore higher absolute rapidity. It should also be noted that for both the 
longitudinal momenta and rapidity predictions, the ME correction method underestimates 
the amount of radiation in the tails of the distributions in comparison to the MCONLO plots. 
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5. Conclusions 



In this work we have realized the MC@NL0 matching prescription for the decay of Higgs 
bosons produced in association with vector bosons at both the ILC and hadron collid- 
ers. This work was achieved within the framework of the Herwig++ Monte Carlo event 
generator. 

We compared the MC@NL0 predictions with those obtained via the matrix element 
correction method as well as leading order predictions. The effects of the hard radiation 
are visible in the reconstruction plots for the mass and energy of the bb pairs resulting from 
Higgs decay. Also visible is the effect on the longitudinal momenta and rapidity of the bb 
pairs. 

Less visible is the effect on the pr spectra for qq' annihilation due to the dominant 
effects of initial state radiation. For e + e~ annihilation, the MCSNLO method predicts a softer 
spectrum for the pt with respect to the beam axis and a harder spectrum with respect to 
the associated Z boson. 

The algorithm used for these processes will be publicly available with the forthcoming 
version of Herwig++. Although we have considered the associated production of Higgs 
and vector bosons in this paper, the MC@NL0 algorithm can be interfaced with other Higgs 
boson production mechanisms, since here we only apply the method to the decay process. 
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A. Monte Carlo algorithm 



The integrals in (3.12) can be evaluated using a variety of Monte Carlo methods. In this 
paper, the 'Hit or Miss' Monte Carlo method is used. This is the simplest and oldest form 
of Monte Carlo integration and essentially involves finding the area of a region in phase 
space by integrating over a larger region, a binary function which is 1 in the region and 
elsewhere. The sampling method used for the points xg, xg is the importance sampling 
method whereby more samples are taken from regions where the integrand is large and 
less from regions where it is small. This ensures that the sampled points have the same 
distribution as the integrand. 

The following algorithm summarizes how the starting bb and bbg configurations were 
generated according to equation 3.12. In the discussion that follows, we work in the rest 
frame of a Higgs boson of mass Mh = 114 GeV and as = as(Mn) = 0.114. 



1. Randomly sample points x&,xg, in each of regions Jb, Jg, D and O of the phase 
space and using the 'Hit Or Miss' Monte Carlo method, evaluate the 5 integrals, 



' 2 ^ , I^p , iff , iff and iff as well as their absolute sum, I. 
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Note also the maximum values of the integrands in I^p and iff. 



r(3)i I/O) I 

and — r — respectively. 



2. The eventual proportion of bb Monte Carlo events will be determined by the ratio 

|/ (2) | + |/ (2) | + |/ (2) | 

J j — . Likewise, the proportion of bbg events in the soft regions Jb, Jb 

and the hard region D are determined by the ratios — — o.m — 
The algorithm below is then used to importance-sample the bbg events so that the 
corresponding xg) values of the Monte Carlo events have the same distribution 
as the integrands in iff andl^ : 

(a) For event generation in region L (where L is one of D, Jb or Jg), randomly 
select a point xg, xg in that region. 



(3) 

(b) Evaluate the absolute value of the integrand in I L for this point, | w{x\ )J xj ) 



Is | w(xb,xi) I > R 



Wv, 



(R is a random number between and 1 and 
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I u>ma,x | is the maximum value of | w(xb,xi) | determined in Step 1). 

(c) If NO, return to (a). If YES, accept the event and set w unw = sgn w(xf J ,xi) i.e. 
w unw = 1 if w(xb,xi) is positive and —1 if negative. (In regions Jb and Jg, 
M < Mc and hence the integrands and the integral, if^ in these regions are 
negative). This process is called unweighting. 

(d) Repeat the process until the correct proportion of bb and bbg events have been 
generated. 

(2 3) 

(e) Using the importance-sampled points, obtain an estimate for the integral, I L = 

— x /, where N is the total number of Monte Carlo events generated. We 
typically use N = 10 6 . 



This is the MC@NL0 method. In this way, for a total of N events, the correct proportion 
of bb and bbg events with ± unit weight is generated with the same distribution as the 
integrands in (A.l). All of these integrals are finite, but the integrands are divergent at 
isolated points within the integration regions. Before the sampling could be carried out, 
the divergences in the integrands (which cause problems in the sampling process) had to 
be taken care of. This is the described in section B. 



B. Divergences and mappings 



B.l Divergences in dead region 

In region D, the hard matrix element squared M given in equation 2.11, diverges as 
(xb,xi) — > (1, 1). To avoid this divergence, one can map the divergent region into another 
region in such a way that the divergence is regularized. This is ensured by the fact that the 
region of integration vanishes as the singularity is approached. There is a double pole in M. 
at (xb, xi) = (1, 1). To avoid this pole, the region Xb, > | is mapped into a region which 
includes D but whose width vanishes quadratically as Xb,x b — > 1 [8,23]. The mapping is: 



x b = 1 - 



\-{l-x b ) 
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- A -x b . 



Xt 



l-2(l-x b ) 



4 V 



%b) 



5 1 3 

8 + 2 Xb+ 2 X ~ b 



2xbXh 



(B.l) 



in the 



when Xb > x b > |. This mapping also introduces an extra weight factor of 2(1 — 
integrand. Interchange Xb and x b in both the mapping and weight factor when x b > Xb > §. 
Figure 14 shows the region mapped (solid) and the region mapped onto (dashed). 

B.2 Divergences in jet regions J& and J5 

In both regions Jb and Jg, there is a simple pole in the term (A4 — Mc) a t (xb, £5) = (1, 1). 
In the region Xb,x b > |, a new set of random points are generated which have a weight 
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Figure 14: The mapped region (solid) and the region mapped onto (dashed) 

factor to cancel the divergence. The mapping used in region Jb where xg > Xb is [8]: 

x' b = l- 0.25ri , 

x'- b = l-(l- x ' b )r 2 (B.2) 

where r\ and r 2 are random numbers in the range [0, 1]. The weight factor for this mapping 
is 27"i . For region Jg, where Xb > xg, interchange Xb and xg in the mapping. The mapped 
regions are shown with solid boundaries in Figure 15. 
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Figure 15: Mapped regions 
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